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SHAKEN AND STIRRED: CONDUCTION AND TURBULENCE IN CLUSTERS OF GALAXIES 



0^ 

o 
o 

o 

Q 



O 

U 

6 

(N 
> 
oo 

in 



o 



M. RUSZKOWSKI^ 

Department of Astronomy, University of Michigan, 500 Church Street, Ann Arbor, MI 48109, USA; 

e-mail: mateuszrOumich.edu (MR) and 
The Michigan Center for Theoretical Physics, 3444 Randall Lab, 450 Church St, Ann Arbor, MI 48109, USA 

S. Peng Oh ^ 

Department of Physics, University of California, Santa Barbara, CA 93106, USA; 
e-mail: peng@physics.ucsb.edu (SPO) 
Draft version December 16, 2009 

ABSTRACT 

Uninhibited radiative cooling in clusters of galaxies would lead to excessive mass accretion rates 
contrary to observations. One of the key proposals to offset radiative energy losses is thermal 
conduction from outer, hotter layers of cool core clusters to their centers. However, thermal 
conduction is sensitive to magnetic field topology. In cool-core clusters where temperature decreases 
inwards, the heat buoyancy instability (HBI) leads to magnetic fields ordered preferentially in the 
direction perpendicular to that of gravity, which significantly reduces the level of conduction below 
the classical Spitzer-Braginskii value. However, the cluster cool cores are rarely in perfect hydrostatic 
equilibrium. Sloshing motions due to minor mergers and stirring motions induced by cluster galaxies 
or active galactic nuclei (AGN) can significantly perturb the gas. The turbulent cascade can then 
affect the topology of the magnetic field and the effective level of thermal conduction. We perform 
three-dimensional adaptive mesh refinement magnetohydrodynamical (MHD) simulations of the effect 
of turbulence on the properties of the anisotropic thermal conduction in cool core clusters. We show 
that very weak subsonic motions, well within observational constraints, can randomize the magnetic 
field and significantly boost effective thermal conduction beyond the saturated values expected in 
the pure unperturbed HBI case. We find that the turbulent motions can essentially restore the 
conductive heat fiow to the cool core to level comparable to the theoretical maximum of '--^ 1/3 
Spitzer for a highly tangled field. Runs with radiative cooling show that the cooling catastrophe can 
be averted and the cluster core stabilized; however, this conclusion may depend on the central gas 
density. Above a critical Froude number, these same turbulent motions also eliminate the tangential 
bias in the velocity and magnetic field that is otherwise induced by the trapped y-modes. Our results 
can be tested with future radio polarization measurements, and have implications for efficient metal 
dispersal in clusters. 
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1. INTRODUCTION 

X-ray clusters show a strong central surface bright- 
ness peak, and many of them have short central cooling 
times. However, Chandra and XMM-Newton observa- 
tions have shown that actual gas cooling rates fall signifi- 
cantly below classical values from a standard cooling fiow 
model; emission lines such as Fe XVI I expected from low 
temp erature gas are not observed (jPeterson fc Fabianl 
l2006f ). This suggests that besides quasi-hydrostatic equi- 
librium, clusters might also be in quasi-thermal equi- 
librium, with some form of heating balancing cooling. 
Many forms of heating have been postulated. Hybrid 
models can account for the observations; leading con- 
tenders generally involve some combination of AGN heat- 
ing and thermal conduction. The remarkable discovery 
by Chandra of AGN blown bubbles in clusters has trig- 

^ work performed in part while on leave at Max Planck Institute 
for Astrophysics, Karl-Schwarzschild-Str. 1, 85741 Garching, Ger- 
many and at Institute of Astronomy, Madingley Road, Cambridge 
CBS OHA, United Kingdom 



gered a widespread renaissance in the study of AGN 
feedback (see MpN amara & Nulsen (2007) for a recent 
review). It has been suggested that the AGN feedback 
alone, generally fails to distribute heat in the spatially 
distributed fashion required to offset coohng (e.g., Ver- 
naleo & Reynolds 2006), both in its radial dependence 
and solid angle coverage. However, recent simulations by 
Briiggen & Scannapieco (2009) involving subgrid model 
for Rayleigh- Taylor-driven turbulence show that AGN 
heating can self-regulate cool cores. In these non-MHD 
simulations the subgrid turbulence model is crucial to 
simulate the evolution of the cool core. In a realistic sit- 
uation, the ability to self-regulate undoubtely depends on 
a range of the mechanisms used to transmit the energy 
from the AGN blown bubbles to the ICM (e.g., dissipa- 
tive sound waves, pdV work, cosmic rays, suppression 
of both the Rayleigh- Taylor instability and the mixing 
of the bubble material by (even weak) magnetic fields). 
These issues are a matter of intense debate. At the same 
time, it is a remarkable fact that if one simply uses ob- 
served temperature profiles in clusters to construct the 
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Spitzer conductive flux from the cluster outskirts, it is 
very nearly equal to that required to balance the ra- 
diative cooling rate as indicated by the observed X-ray 
surface brightness profile, for some reasonable fraction 
/ ~ 0.3 of the Spitzer value (e.g.. Fig. 17 of Peterson & 
Fabian 2006). There is no reason in principle why such 
close agreement should exist, and it is a tantalizing hint 
that nature somehow "knows" about Spitzer conductiv- 
ity. Models which invoke both mechanisms are able to re- 
produce observed temperature and density profiles with- 
out any fine-tuning of feedback or conduction param- 
eters (Ruszkowski & Bcgclman 2002; Guo & Oh 2008), 
unlike conduction-onlv models (jBertschinger fc MeiksinI 
[19861: iZakamska fc NaravMl[2003[ r 

Nonetheless, even if no fine-tuning of conductivity is 
required, it is unsatisfying to invoke an ad-hoc suppres- 
sion factor /; the anisotropic transport of heat along 
field lines can be calculated from first principles. If the 
field lines are t angled and chaotic, then / ^ 0.2 might 
be appropriate (jNaravan fc Medvedevll200"H ). However, 
recent work has uncovered instabilities arising in strat- 
ified plasmas, which rearrange the field lines on large 
scales, with strong implications for heat transport. Un- 
like convective instability, which only arises if entropy 
declines with radius, these instabilities depend on tem- 
perature gradients, and arise irrespective of the sign 
of the gradient. If dT/dr < 0, as in the outer re- 
gions of clusters, the resu lting magnetothe rmal instabil- 
itv (M TI; Balbus (200(1): IParrish fc Stond (|2005. 2007); 
[Parrish ct al. (2008)) drives field lines to become pref- 
erentially radial, substantially boosting conductivity to 
close to the Spitzer value, / ~ 1. On the other hand, 
if dT/dr > 0, as in the cores of cool-core clus ters, then 
the re s ulting heat bu o yancy instability (HBI; lOuataertI 
(pool : IParrish et all (j2009l ). Bogdanovic et al. 2009) 
will drive field lines to become preferentially tangential, 
essentially shutting off radial conduction (/ ^ 1), and 
leading to catastrophic coohng in the cluster core. Thus, 
it would seem that a more careful calculation vitiates the 
ad-hoc assumptions of models which invoke conduction. 

However, this presumes that there are no compet- 
ing mechanisms which also affect field line topology. 
One possibilit y which has been speculated upon in 
the literature (iRuszkowski et all [20071: iGuo et al.l [20081 : 
iBogdanovic et all I2009D is AGN bubbles, which upon 
rising would leave radia lly orientated field lines in their 
wake. iGuo fc OhI ([2009') showed that time-variable con- 
duction, perhaps triggered by AGN outbursts, could al- 
low clusters to cycle between the cool core and non- 
cool core states. Another, perhaps more general mech- 
anism is simply turbu lent gas motions, which ca n 
be induced by mergers (Ascasibar fc MarkevitcW 12006*) , 
the m otion of galaxies (Kim 2007: Conrov fc Ostrikcr 
[2008'), AGN outbursts (McNamara & Nulscn 2007), or 
cosmic-ray driven convection ( Chandran fc Rasera,2007l : 
iSharma et gd] l2009a[ )^. Indeed, iSharma et gd] (|2009bl ) 
suggested that the amount of turbulence needed to avoid 
the stable end-state induced by the HB I could bees- 
timated from the Richardson number ([Turneil Il973f ). 

^ Note, however, that since the HBI is a buoyancy instability 
rather than a convective instabihty, turbulence cannot be self- 
consistently generated by the HB I itself, since the induced veloci- 
ties are small (.Bogdanovic et al.ir2009i : , Parrish et al.,,2009.) . 



which is simply the ratio of the restoring buoyant force 
to the inertial (pu • Vu) term. If one defines Ri = 
gr{d\nT / d\nr) / for typ ical cluster conditions , the crit- 
ical turbulent velocity is (jSharma "eraI]l2009bD : 
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where g-s is the gravitational acceleration in units of 
10~^cm^s~^, rio is a characteristic scale height in units 
of 10 kpc, and Ric is the critical Richardson number. 
Ric ^ 1/4 is typical for hydrodynamic flow; order unity 
differences might exist for the MHD case. Such lev- 
els of turbulence are easily seen in cosmologica l sim- 
ulations, even in very relaxed clusters ( Evrard 



Norman fc BrvanI 119991: iNagai et all l2003t iVazza et all 
2Q09a|[b!), presumably due to continuing gas accretion 
and the supersonic motions of galaxies through the ICM. 
Gas "sloshing" in cl usters has been invoked to explain ob- 
served cold fronts ([Ascasibar fc MarkevitchI [20061 . and 
biases in X-ray ma ss profiles, particularly in the clus- 
ter core (|Lau et al . 2009). In the future, direct mea- 
surement of turbulent velocities could best be performed 
by a high-spectral resol ution imaging X-ray spectrome- 
ter ([Rebusco et al.l[2008 ), although early upp er limits ex- 
ist from XMM-Newton ([Sanders et al.l ((20090 : for A1835, 
non-thermal broadening is less than ~ 275kms~^). Be- 
sides the strong impact on conductive heat transport, the 
interaction of turbulence with conduction has testable 
implicati ons for metal disper sal (which is strongly en- 
hanced; ([Sharma "etalll2009ari '). With an instrument as 
as the Square Kilometer Array, the field topology itself 
can potentially be mapped out via the rotation measure 
of background radio sources. 

Our goal in this paper is to quantitatively examine for 
the first time if indeed turbulent motions can overwhelm 
the tangential field topology driven by the HBI (we defer 
the MTI dominated regime to future work), and allow 
thermal conduction at the level required to avert a cool- 
ing catastrophe. We begin with a MHD FLASH simu- 
lation of a cluster which in its unperturbed state suffers 
the HBI, consistent with previous work. We then simu- 
late turbulence via a spectral forcing scheme that enables 
statistically stationary velocity fields (Eswaran & Pope 
1998). An important feature of our study is that even if 
the medium is stirred isotropically, the turbulent velocity 
field itself can be anisotropic, and interact in non-trivial 
ways with thermal conduction. This is because when the 
turbulent driving frequency falls below the Brunt- Vaisala 
frequency, (^modes can be excited; the trapped modes 
induce tangentially biased vortices, due to the r estor- 
ing forces which act on vertical motion (Rilcv & Lelongj 
[2000.') . Thus, even in the absence of the HBI, the excita- 
tion of (/-modes can tangentially bias the magnetic field. 
In realistic clusters, (;-modes also greatly enhance the 
efficacy of stirring, since trapped modes allow the tur- 
bulence to be volume-filling, although our present simu- 
lations already have volume-filling turbulence by design. 
The Brunt- Vaisala frequency, which determines where 
g-modes can be produced, depends both on the gravi- 
tational potential and the gas temperature and density 
profile. The latter can in turn be affected by heating 
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and thermal conduction. Thus, the cluster potential, gas 
properties, and thermal conduction interact in often sub- 
tle ways. Overall, we find that turbulence can indeed 
restore conductive heat flow back to the values / ~ 0.3 
expected for a highly tangled field, and required to stem 
a cooling flow. 

The organization of this paper is as follows. In Sj^l we 
describe simple scaling relations to guide our intuition; 
in fj3l we describe our computational methods; in 21 our 
results; we then conclude in fj5l 

2. THEORETICAL EXPECTATIONS: TURBULENCE IN A 
STRATIFIED MEDIUM 

It is useful to begin by considering theoretical expec- 
tations for fluid motions in a stably stratified medium. 
We first review the case for a hydrodynamic fluid, which 
in clusters is stabilized by a strong entropy gradient. We 
then consider the MHD case where anisotropic conduc- 
tion alters the picture. 

In the purely hydrodynamic case, the medium is con- 
vectively unstable if d\nS/d\nr < 0. On the other hand, 
if the medium is stably stratified with d\nS/d\nr > 0, 
a fluid element which is adiabatically displaced about 
its equilibrium position will experience a restoring force 
f ~ pr = —pg{3/5){dhiS/dr)Ar, assuming it remains 
in pressure balance with its surroundings. This is the 
equation for simple harmonic motion, at frequency: 



^BV 



g 3 d InS 
r 5 dlnr ' 



(2) 



the Brunt- Vaisala frequency. A formal WKB analysis re- 
veals that the dispersion relation for such ji-modes, where 
gravity is the restoring force, is: 



-'BV 



(3) 



where fc^ = k'^ + kf.. This immediately then implies that 
for oscillatory modes with real kr, uj < i^BV] otherwise 
the modes have imaginary radial wavenumber and are 
evanescent. Physically, one can always achieve a low 
frequency response by making the mode progressively 
more tangential, but it is impossible to drive the sys- 
tem at frequencies higher than the maximum response 
frequency of wbvj corresponding to completely vertical 
oscillations. This thus implies that waves driven at fre- 
quencies UJ < uj-BV can be resonantly excited, and will 
be trapped, reflected and focused inside the resonance 
radius where Ci^BV = w. This has implications for how 
stirring can be amplified into volume-filling turbulence. 

For our purposes, an even more interesting aspect 
of such g-modes is that they tend to be preferentially 
tangential, as has been widely confi rmed in calcula- 
tions of stellar pu lsation (jCoxi llQSOD . galaxy clusters 
(ILufkin et al.l 119951). and th e earth's atmosphere and 
oceans ( Rilev fc Le long 2000) . Let us consider a simple 
order of magnitude argument for why this should be the 
case. For simplicity, consider an isothermal atmosphere. 
In steady state, and adopting the Boussinesq approxi- 
mation (such that V • V = 0), the continuity equation 
reads v • Vp = 0. Since the background density po does 



not vary horizontally, the perturbed continuity equation 
is V • Vp' +Vzdpo/dz = 0. If |v| ~ U, \vz\ ~ W, and L is 
a characteristic lengthscale, then: 



WL 
IT 



dpo 
dz 



(4) 



The buoyancy forces from stratification impose addi- 
tional constraints. The vorticity evolution form of 
the momentum equation'^ for waves (i.e., assuming 
SP/P < 6p/p) is (Lufltin ct al. 1995.): 



dt 



= i — {k X g) 
P 



(5) 



which gives to order of magnitude, 



P ^ Po—r- 

Combining equation ([5]) and ([7]), we obtain: 



(6) 



W 
If 



PoU^ 



( u 



gL'^\dp/dz\ \lubvL 



Fr^ 



(7) 



where Fr is the Froudc number, which in this case com- 
pares the characteristic frequencies of inertial and grav- 
itational (buoyancy) forces. Ri ~ 1/Fr'^ is often known 
as the Richardson number. Equation ([5]) reveals two 
important properties. When turbulence is absent or 
very weak, Fr <ti 1 , motion is fundamentally horizontal, 
W <^ U: strong restoring buoyancy forces oppose mo- 
tion in the vertical direction. However, stronger stirring, 
Fr <; 0{1), can overwhelm such forces and re-establish 
isotropic turbulence. 

While all of the preceding discussion hold only for short 
wavelength modes in the WKB approximation, numeri- 
cal simulations of non-linear evolution of global g-modes 
in galaxy clusters excited by orbiting galaxies show that 
they can still be very usefully interpreted in terms of the 
linear WKB analysis (jLufkin et al.lll995l ). 

How do these expectations change in the presence of 
anisotropic conduction? Here, behavior is governed by 
temperature rather than entropy gradients, and the at- 
mosphere is buoyantly unstable regardless of the sign 
of VT. For the HBI, the instability saturates once the 
magnetic field is perpendicular to VT. Displaced fiuid 
elements now experience a restoring force f = pr = 
—pg{d\nT/dr)Ar (Sharma et al. 2009b), and hence os- 
cillate about their equilibrium position at a frequency: 



,-, ,mhdn2 _ 9 dlriT 
l^BV ) - --7j — ) 
r dlnr 



(8) 
(9) 

which is typically ~2 times smaller than ujbv- Note that 
this assumes tcond ~ L'^/k <C isc ^ L/cs (where k ~ 

^ Note that since g-modes generically generate vorticity, vorticity 
maps are a very useful marker of their presence. 
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vJe is the conductivity coefRcient) so that a displaced 
blob's temperature is determined by conductive rather 
than adiabatic cooling. This requires: 



200 kpc 



tkcV j 



(10) 



which in any case is required by t he WKB assum ption, 
L <C r. The dispersion relation is (|Quataertil2008[ ): 



/ MHD\2 



(1-26^)1^ + ^^4^^ , (11) 



fc2 



(12) 

where z = f, bx,z = Bx.z/B and B = S^^x + BzZ, 
without loss of generality. One can then solve the 
quadratic equation for kz to obtain the criterion for os- 
cillatory solutions (real kz), which yields the criterion 



w ^ (1 - 2h 



2\l/2, ,MHD 
' '^BV 



,MHD 
■^BV 



(1/ \/3, 1) , correspond- 



ing to fully tangled (tangential) fields. Since temper- 
ature gradients are generally milder than entropy gra- 
dients, the criteria for g-modes is therefore marginally 
more stringent in the MHD case. The discussion of tan- 
gentially biased motion carries over with equal force to 
the MHD case, except for the change in the nature of the 
restoring force, such that: 



W 
U 



U 



,MHD 
"^BV 



L 



(13) 



so that if Fr^^^ <; 0(1), turbulent motions can 
isotropize fluid motions (and thus the magnetic field). 

We now proceed to exphcitly test these notions in 3D, 
MHD simulations. It is important to note that uj^y^ 
itself can be appreciably altered by gas motions. For 
instance, if thermal conduction is restored by stirring, 
heat transfer acts to reduce the temperature gradient 
and hence w^y^, further increasing Fr^^^ _ In the 
limit where efficient conduction causes the cluster to 



become isothermal, 



,MHD 
^BV 



and fluid elements 



are neutrally buoyant, with no restoring force upon 
displacement. Such a situation may indeed prevail in 
clusters which have fairly flat temperature profiles. In 
cool core clusters, the temperature profile likely depends 
on a complex interplay between turbulence, conduction, 
radiative cooling and AGN feedback. 

3. METHODS 
3.1. Initial conditions for the gas 

For concreteness, we adopt simulation parameters ap- 
propriate for the well-observed Tvir ''^ 5 keV cool-core 
cluster A2199. Note that although we focus attention on 
the cluster core in the paper, we simulate the entire clus- 
ter in a large 1 Mpc box, to avoid sensitivity to assumed 
boundary conditions. 

We model the dark matter distribution as a softened 
NFW profile of the form: 



Pdm — 



Mo/27: 



{r + rc){r + 



(14) 
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Fig. 1. — Initial density (solid line) and temperature (dashed 
line) in the ICM. 



where — 20 kpc is the smoothing core radius and = 
390 kpc is the usual NFW scale radius. The parameter 
Mo = 3.8 X IO^'^Mq determines the cluster mass and 
is of the order of the cluster virial value. This density 
distribution corresponds to a gravitational acceleration: 



(r - rcY 



rs{rs - r-c) 



r,(r, - 2r,) In ( 1 + - ) + In ( 1 + - ) (15) 



To model the gas, we adopted a power law entropy pro- 
file w ith a core, as seen in observations (jCavagnolo et al.l 
[2009h : 



K{r) ^Ko + Ki 



i- 

\r200 



(16) 



where entropy K = kBT/m?^^. Here we adopt a — 1.1 
and adjust Kq and Ki to match observed temperature 

and density profiles for A2199 {Kg = ksTo/riQ^^ where 
Tq, uq are observed central temperature and density, and 

2/3 

Ki — K20Q = fcBT2oo/n2oo; we estimate from the virial 
theorem fcsT2oo « GM2oo'Tjp/2r2oo), and assume that 
the gas profile traces the dark matter profile in the outer 
regions. Given this entropy profile, we can solve for hy- 
drostatic equilibrium: 



dlogp 
d logr 



d logs' GM 



dlogr rkgT 



(17) 



where ksT — Sn"^^^, to obtain the temperature and den- 
sity profiles. These initial conditions can be seen in Fig- 
ure [1 
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3.2. Magnetic fields 

We generated stochastic magnetic fields using a 
method similar to that described in Ruszkowski et al. 
(2007). In brief, we set up magnetic fields in Fourier 
space and generate independent field components in all 
three directions in k— space. This is done to ensure that 
the the field has a random phase. For example, for the 
a:-component, we set up a complex field such that 

(Re(B,(k)),Im(B,(k))) = {G{ui)B,Giu2)B), (18) 

where G is a function of a uniformly distributed indepen- 
dent random variables ui or U2 that returns Gaussian- 
distributed random values. The amplitude B is given 

by 



B (X fc""/^ exp 



(19) 



where = 27r/Ao and Aq ~ 43/i ^ kpc is a smoothing 
wavelength. We then apply a divergence cleaning opera- 
tor in k— space as follows: 



B(k) — > (l-kk)B(k), 



(20) 



where k is the unit vector in k-space. Note that this 
operator does not change the shape of the power spec- 
trum of the magnetic field fluctuations. This is sufficient 
for our purposes as we assume the field is dynamically 
unimportant. We then use three-dimensional inverse fast 
Fourier transform to convert the field in k-space to the 
real space. 

3.3. MHD equations with anisotropic thermal 
conduction and radiative cooling 

We solve the MHD equations, including field-aligned 
thermal conduction transport. These equations read: 



dp 

at 



V • (pv) = 



g(pv) 
dt 



dE 



V • (pvv - BB) + Vp = /9g 



— +V-(v(£;+p)-B(vB)) = 
pg • V - V • F + C 



5B 



V • (vB - Bv) = 0, 



where 



B^ 



P = Pth 



E ^ i-— +e+ — 



(21) 
(22) 

(23) 
(24) 

(25) 
(26) 



where pth is the gas pressure and e is the gas internal 
energy per unit volume. We assumed the adiabatic index 
7 = 5/3. The anisotropic thermal conduction heat fiux 
F is given by 



F = -KeB(eB- VT), 



(27) 



where is a unit vector pointing in the direction of the 
magnetic field and k is the Spitzer-Braginskii conduction 
coefficient given hy n — 4.6 x 10~^T^/^erg s~^cm~^K~^. 

In equation (21), C represents the cooling rate per unit 
volume. We u se standard tabulated publ icly available 
cooling curves (jSutherland fc DopitalfTOOl for metallic- 
ity Z — 0.3^0. 

Magnetic field evolution was solved by means of a 
directionally unsplit staggered mesh algorithm (USM; 
Lee & Deane 2009). The USM module is based on 
a finite-volume, high-order Godunov scheme combined 
with constrained transport method (CT). This approach 
guarantees divergence-free magnetic field distribution. 
We implemented the anisotropic conduction unit follow- 
ing the approach of Sharma & Hammet (2007). More 
specifically, we applied monotonized central (MC) lim- 
iter to the conductive fluxes. This method ensures that 
anisotropic conduction does not lead to negative temper- 
atures in the presence of steep temperature gradients. 
Tests of the method are presented in accompanying pa- 
per (Ruszkowski et al. 2009, in preparation). 

3.4. Driven turbulence 

In order to emulate the effects of turbulence in the 
ICM, we employed a spectral forcing scheme that en- 
ables statistically stationary velocity fields (Eswaran & 
Pope 1988). This scheme utilizes an Ornstein-Uhlenbeck 
random process, analogous to Brownian motions in a vis- 
cous medium. At each timestep, accelerations are ap- 
plied to the gas. The acceleration field has vanishing 
mean and each acceleration mode has constant disper- 
sion and is time-correlated. That is, at each timestep 
the acceleration in a given direction is a Gaussian ran- 
dom variable with a fixed amplitude and decays with 

time as oc \/l — where / = exp(— t/rdccay)- This es- 
sentially means that the forcing term has a "memory" 
of the previous state. The phases are evolved in Fourier 
space and the divergence in the flow is cleaned making 
the flow incompressible, W ■ v ~ 0. This is consistent 
with the Boussinesq approximation, and justified since 
the gas motions are significantly subsonic. Further de- 
tails and the numerical tests of this method can be found 
in Fisher et al. (2008). The key quantity of interest is 
the velocity dispersion of the resultant velocity field a 
and its scaling with the injected energy and the injection 
spatial scale. This can be understood in terms of sim- 
ple dimensional analysis in an equilibrium situation when 
the turbulent driving is balanced by viscous dissipation 
in the fluid 



(28) 



where iVmodc is the number of Fourier driving modes, 
e* is the energy injection per unit mass per mode and 
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V ~ Lv is the gas viscosity for a a given largest eddy 
scale L. This implies the following scaling of the velocity 
dispersion 

a e(A^modci)'/'e*'/', where (29) 

We tested this scaling in a set of small box simulations 
and use it as a guide for choosing the parameters for 
the cluster simulations. Table 1 provides a summary of 
the turbulence parameters. The values of the velocity 
dispersion quoted there are for the conductive case 
averaged within ^ 100 kpc from the cluster center. For 
closed box simulations, the velocity field is approxima- 
tion Kolmogorov, though this is of course altered in a 
stratified medium. 

The simulations were performed using the FLASH 
code (version 3.2). FLASH is a modular, parallel 
magnetohydrodynamic code that possesses adaptive 
mesh refinement capability. The three-dimensional 
computational domain was approximately 1 Mpc on 
each side, enclosing a large fraction of the cluster. 
The central regions of the cluster had an enhanced 
refinement level. The maximum spatial resolution for 6 
levels of refinement was ~ 2.7h~^ kpc. We performed a 
set of shorter simulations to verify that the results did 
not depend on the resolution for the chosen set of initial 
condition parameters. The simulations were performed 
on a 320-processor cluster located at the Michigan Aca- 
demic Computing Center at the University of Michigan 
in Ann Arbor and on the Columbia supercomputer at 
NASA Ames. 



4. RESULTS 

In Figure 2 we present cross-sections through the 
cluster center. Both panels show the orientation of 
magnetic field vectors in the final state. Superposed 
on these vectors is the gas density distribution. Left 
hand side panel corresponds to the HBI case and the 
right one to the anisotropic conduction in the presence 
of weak stirring gas motions. In the latter case, the 
mean velocity dispersion is significantly subsonic and of 
the order of ^ 50 km/s. It is evident from this figure 
that the addition of weak stirring "outcompetes" HBI 
and dramatically changes the topology of the field. As 
explained below in detail, this randomization affects the 
effective level of thermal conduction. 

In Figure 3 we quantify the effect of HBI and tur- 
bulence on the properties of the intracluster medium. 
The parameters of all runs are summarized in Table 
1 and their meaning is explained below. All curves 
correspond to the quantities averaged within '--^lOO kpc 
from the cluster center. In the left panel we show the 
evolution of the anisotropy in the velocity field. The 
anisotropy is defined a,s /3y — 1 — /{2a'^), where at and 
ar are the velocity dispersion in the tangential direction 
and radial directions, respectively. Thus, the isotropic 
case corresponds to Py — 0, while positive (negative) 
values of are for preferentially radial (tangential) gas 
motions. 

We begin by discussing the pure HBI case (dashed 



lines) in all three panels in Figure 3. In this case, the 
initial tangled magnetic field allows thermal conduction 
at the level of ~ 1/3 of the Spitzer value, and the 
resulting transfer of heat leads to radial flows to adjust 
to a new equilibrium. This rapidly shuts off once 
thermal conduction is suppressed. After the initial 
readjustment, the velocities become progressively tan- 
gential (left panel). This is consistent with theoretical 
expectations for the low Froude number case, as in ^J2j 
strong buoyancy restoring forces confine motion to be 
largely tangential. 

In the middle panel we show the anisotropy in mag- 
netic field Pb — ^ — '^Bt/i'^^Br)^ where aBt and a Br are 
the dispersion in the tangential and radial magnetic field 
components, respectively. This figure clearly illustrates 
that the HBI leads to strongly tangentially-biased 
magnetic fields. Finally, in the right panel we present 
the effective intrinsic thermal conduction coefficient 
of the ICM in terms of the Spitzer-Braginskii value. 
More specifically, the intrinsic effective conduction 
coefficient of the ICM is /k, where / is given by 

/ = (-FV)/(i^r,spitzor), whcrC 



(Fr) (X (T5/2|(eT-eB)(eB-e,)|), (30) 

and 



(-Fr.spitzcr) « (r5/2|eT-e,|). (31) 

The volume averaging in the above expressions is per- 
formed within the central '--^lOO kpc. The proportionality 
constant in Eq. 28 and Eq. 29 are identical. The versors 
ex, gb and point in the direction of the temperature 
gradient, magnetic field vector and the radial direction, 
respectively. Note that the projection of the temperature 
gradient onto the B-field is given by the averaged abso- 
lute value, rather than a straight average. This correctly 
reflects the expectation that for = and a com- 
pletely randomized B-field, / — {{es ■ e^)^) = 1/3. We 
refer to / as the Spitzer fraction below. The effective con- 
duction declines from the value close to the theoretical 
maximum ^ 33% to about 4% of the Spitzer-Braginskii 
value. The initial values of this ratio are somewhat lower 
then the theoretical value. This simply reflects the fact 
that we are sampling a finite volume of the cluster and so 
the actual value of the ratio depends on the random seed 
used to generate the field distribution. The decline in the 
conductive fiux is a consequence of the rearrangement of 
magnetic fields to become preferentially tangential. 

We now compare the HBI reference results to simu- 
lations that include driven turbulence and anisotropic 
conduction. The solid black curve shows results from the 
weak turbulence driving case. The total energy injection 
due to stirring motions is much smaller then the thermal 
energy content of the ICM. The gas velocity dispersion 
in the final state is ~ 50 km/s, i.e., much smaller the 
the sound speed in the ICM. Such weak stirring motions 
are entirely conceivable and are consistent with, for ex- 
ample, repeated minor mergers, major mergers that stir 
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TABLE 1 

Turbulence parameters of the performed runs 





k ■ 

"'mm 


k 

"'max 


'^dccay 


e* 






Name 


(cm^^) 


(cm-1) 


(s) 


(cm^/s- 


') 


(km/s) 


weak 


2.0 X 10~23 


5.1 X IQ-^^ 


3.1x10^5 


8.1 X 10 




~50 


strong 


2.0 X 10-23 


5.1 X IQ-^^ 


S.lxlQis 


3.0 X 10 


-9 


~150 




Fig. 2. — Topology of the magnetic field (vectors) with superimposed gas density distribution (color background). The results are shown 
in the plane intersecting the cluster center. The images are ~1 Mpc on a side. Left panel corresponds to the end state of the simulation with 
anisotropic conduction and the right panel to the simulation involving both the anisotropic conduction and very weak subsonic turbulence. 



the gas on large scales and cascade to smaller scales, 
AGN-dri ven motions or stir ring by galaxy motions. For 
instance. INagai et al.l (|2003f ) find from cosmological sim- 
ulations that even relatively relaxed clusters have inter- 
nal velocities ~ 20% of the internal sound speed; iKimI 
dlOO?) finds turbulent velocities 100 - 200 km s'^ in 
numerical simulations of gravitational wakes of galaxies 
in clusters. The left panel demonstrates that the gas ve- 
locity field becomes preferentially tangential in this case, 
while the middle panel shows that the magnetic field is 
close to isotropic. This magnetic field topology is consis- 
tent with the right panel in Figure 2. The isotropization 
of the magnetic field leads to a significant boost in the 
level of the effective thermal conduction as shown in the 
right panel of Figure 3. 

In order to better understand the observed trends in 
the velocity, magnetic field and the effective conduction, 
we performed MHD runs with stirring but without con- 
duction. An example of the result from one of these 
runs is shown as a solid purple line in all three panels in 
Figure 3. It is evident that the stirring leads to preferen- 
tially tangential velocity field. This behavior is very sim- 
ilar to the case when both the stirring and HBI operate. 
The reason for the preferentially tangential gas motions 
in the absence of conduction is that the characteristic 



driving frequency is lower then the local Brunt- Vaisala 
frequency lobv- In this case, gas motions excite gravity 
waves ((/-modes) that become trapped within the radius 
where "turbulence frequency" ~ a/X < lubv, where a 
and A are the gas velocity dispersion and the character- 
istic eddy size, respectively. Consequently, the distribu- 
tion of the orientation of magnetic field versors is also 
tangential (see middle panel). One subtle point worth 
emphasizing is that the magnetic fields have "memory" 
of the the fluid displacement history. In other words, 
for weak fields, magnetic fields behave as "tracer parti- 
cles" and, thus, the middle panel essentially quantifies 
the anisotropy in integrated fluid displacements. It is 
interesting to note that the magnetic field for weak stir- 
ring is close to isotropic in the conductive case, while it 
remains preferentially tangential in the absence of con- 
duction. We argue that this effect can be qualitatively 
understood in terms of the magnitude of the buoyant 
restoring force. The magnitude of the restoring force in 
the magnetized medium with anisotropic conduction de- 
pends on VT whereas for the pure MHD the restoring 
force depends on VS*, where S is the gas entropy. The 
entropy gradient is steeper than the temperature gra- 
dient, so the buoyant restoring force in the conductive 
case is weaker then in the pure MHD case. This means 
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that it is easier to perturb and isotropizc^ such a fluid. 
Tliis may explain why the distribution of magnetic field 
in a turbulent magnetized and conducting ICM is more 
isotropic than in the absence of conduction. We also 
note that the tangential bias has been seen in non-MHD 
adiabatic cosmological simulations that did not involve 
thermal conduction (Rasia et al. 2004). In these simu- 
lations the characteristic velocity dispersion in the ICM 
was also highly subsonic and the volume-averaged veloc- 
ity field anisotropy (coincidentally defined the same way 
as here) was — 2 < /3 < 0. This result is consistent with 
the idea of trapping of internal gravity waves. 

The results for the strong stirring case are presented 
as solid blue (turbulence without thermal conduction) 
and solid green (turbulence with thermal conduction). 
In these cases, the gas velocity dispersion is of the or- 
der of cr ~ 150 km/s (see next paragraph for the dis- 
cussion of the velocity dispersion). Such motions, while 
being clearly subsonic, are nevertheless sufficiently vig- 
orous to isotropize the velocity field with and without 
anisotropic conduction. They are also sufficiently pow- 
erful to overwhelm the buoyant restoring force and make 
the magnetic field distribution isotropic in both cases. 
The Spitzer fraction is effectively 1/3 in this case (note 
that in the purely hydrodynamic case, the effective con- 
duction is computed at the post- processing stage). 

In Figure 4 we show the evolution of the volume- 
averaged velocity dispersions (left panel) and the profiles 
of velocity dispersions in the final state (right panel). As 
expected, the unperturbed HBI case leads to very small 
velocity dispersions whereas the volume-averaged veloc- 
ity dispersion saturates at ^ 50 km/s and ~ 150 km/s 
in the weak and strong stirring cases, respectively. The 
velocity dispersion profiles also show that the typical ve- 
locities are always significantly subsonic throughout the 
ICM. Black and green lines are for the conductive case 
and purple and blue for the non-conductive ones. 

In Figure 5 we present characteristic frequencies in the 
cluster: orbital (black) , Brunt- Vaisala (blue) wbv for hy- 
drodynamic case, ° (green) for MHD case, stirring 
frequency in the non-conductive case (yellow), stirring 
frequency in the MHD case (red) . The stirring frequency 
is defined as ^ o"/A, where A = 20 kpc is some rcfcr- 
enc;e length of the order of the coherence scale of the 
turbulent velocity field. Left panel is for the weak turbu- 
lence and the right panel for the strong turbulence case. 
This figure is helpful in explaining the trends seen in the 
topology of the field and the Spitzer fraction (see Figure 
3). Specifically, it can be seen that for the weak stirring 
and non-conductive case, the Brunt- Vaisala frequency 
exceeds the stirring frequency. This implies trapping of 
internal gravity waves and tangential bias in the velocity 
and magnetic fields. This is consistent with the mag- 
netic field topology in this case (purple curve, middle 
panel in Figure 3) that shows preferentially tangential 
magnetic fields. In the weak stirring conductive case, 
the stirring frequency is either comparable to the "mag- 
netic Brunt- Vaisala frequency" or it exceeds this 
frequency. This implies that it should be easier to ran- 
domize the magnetic fields. This is also consistent what 
is seen in Figure 3 (black curve, middle panel) that shows 
more isotropic magnetic fields in this case. The results 
for the corresponding strong stirring case are shown in 



the right panel of Figure 5. In both the conductive and 
non-conductive cases, the stirring frequencies exceed wbv 
and iiJ^y^ which, consistent with our findings, leads to 
isotropic magnetic field distribution. 

We now discuss the runs with radiative cooling. We 
performed a number of simulations that simultaneously 
include anisotropic thermal conduction, radiative cooling 
and different levels of stirring. These runs were done for a 
range of gas densities and temperatures. We now present 
two cases for the initial conditions described in Section 
3 and then briefiy comment on the cooling runs for dif- 
ferent background ICM profiles. As above, we consider 
exact same numerical setup as for the previous adiabatic 
cases but now add radiative cooling. As before, the weak 
stirring case has an extremely small turbulent gas veloc- 
ity dispersion a ^ 50kms^^, which is much smaller than 
the gas sound speed. In general, the presence of radiative 
cooling should not dramatically change the conclusions 
related to the topology of the magnetic field. The cooling 
catastrophe is delayed with respect to the initial cooling 
time by a factor of ~ 4.3. Note that the energy in tur- 
bulent motions is much less than the thermal energy of 
the gas. This ratio is: 

q^(^)\lM^ (32) 

where M. is the Mach number defined here as = a/cs, 
and we have used = ^P/p = 'y<7^as/3. For example, 
for T = 3 keV and fj = 50 km/s (cf^. Figure 4), we get 
g ~ 1.4 X 10^"^. From Figure 4, we see that the turbulent 
energy is injected on a timescale 7 x 10^ years, which is 
comparable to the central cooling time. Since g -C 1, 
the turbulent energy injection rate is much smaller then 
the cooling rate. This estimate is independent of the 
presence of radiative cooling. The results for the strong 
stirring equivalent of the above run {q ~ 1.3 x 10^^) show 
that the cooling catastrophe can be averted. The tem- 
perature profile flattens smoothly with time as a result 
of thermal conduction and mixing (see Figure 6).* 

In simulations involving radiative cooling and turbu- 
lence resulting from acoustic-gravity waves, Fujita et al. 
(2004) found that cooling catastrophe can be prevented 
even in the absence of thermal conduction. However, the 
level of turbulence in those simulations was higher then 
in the cases that we consider. 

Figure 7 is the analog of Figure 3 (right panel) and 
shows the behavior of the Spitzer fraction (solid line). 
As in the adiabatic case, the ratio evolves asymptoti- 
cally toward an approximately constant level but in this 
case the saturated level is slightly lower (~ 30% of the 
Spitzer- Braginskii value). The unperturbed HBI case is 
shown for reference (dashed line). We attribute the small 
offset between the theoretical maximum Spitzer fraction 
of 1/3 and the simulated one to the fact that the HBI is 
marginally more competitive in the cooling case, due to 
the steeper temperature profile. 

* We point out that independent simulations involving radia- 
tive cooling were performed by Parrish et al. (2009) who obtained 
qualitatively similar conclusions. In particular, they find that suf- 
ficiently strong stirring can prevent excessive cooling. Our strong 
stirring case with cooling also confirms that finding. 
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1.0x10'' 1.5x10'' 
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Fig. 3. — The evolution of velocity field anisotropy (left panel), magnetic field anisotropy (middle panel) and the effective thermal 
conduction as a fraction of the mean Spitzer-Braginskii conduction, all averaged within the central 100 kpc. Dashed lines correspond to the 
unbridled HBI instability. Solid black and purple lines are for weak turbulence with and without thermal conduction, respectively. Green 
and blue lines correspond to stronger turbulence with and without thermal conduction, respectively. See text for more details. 
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Fig. 4. — Left panel: Evolution of the gas velocity dispersion in the central 100 kpc with time. Right panel: radial profile of velocity 
dispersion. The color coding of the lines is the same as in Figure 3. 



We note that the cooling results described above come 
with possible caveats. The above examples of cooling 
runs may not carry over to higher ICM densities. We 
performed a limited number of numerical experiments 
for higher ICM densities and found that the develop- 
ment of thermal instability is more likely in those cases. 
If so, it is possible that at higher ICM densities AGN 
may be required for stability. This is consistent with the 
findings of Conroy & Ostriker (2008) who, based on ID 
models, suggest that cooling catastrophe should occur 
at densities around 0.02 cm~'^ (somewhat below our ini- 
tial central density, see Figure 1), if no AGN feedback is 
present. The precise value of the critical density thresh- 
old (if confirmed in three-dimensional simulations) would 
have to be determined by detailed numerical simulations. 
We defer this study to future work. 



5. DISCUSSION 

We have shown that a very low level of turbulent per- 
turbations, '--^ 50 — 150km/s to the ICM — entirely 
consistent with the expectations for cosmological infall, 
galaxy motions, mergers, or AGN activity — can en- 
tirely alter the magnetic field distribution resulting from 
the HBI instability. Instead of preferentially tangential 
magnetic fields, the final field configuration can be easily 
randomized. As expected from simple theoretical con- 
siderations, this happens when the typical Froude num- 
ber (equation [HI) Fr^^^ ^ 0(1). Note that weak, low 
frequency perturbing motions lead to trapped g-modes 
and result in preferentially tangential gas motions, due 
to strong restoring buoyancy forces in the vertical di- 
rection; thus, magnetic fields could in principle become 
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Fig. 5. — Characteristic frequencies in the cluster (in [s ^]): circular orbital (black), Brunt- Vaisala (blue) wbv for hydrodynamic case, 
"^BV^ (green) for MHD case, stirring frequency in the hydrodynamical case (yellow), stirring frequency in the MHD case (red). See text 
for definitions. Loft panel is for the weak turbulence and the right panel for the strong turbulence case. 
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Fig. 6. — Time sequence of temperature profiles for strong stirring 
with both thermal conduction and radiative cooling. The black 
curve denotes the initial state, whilst other curves are plotted every 
~ 1 Gyr. Unlike the HBI case, no cooling catastrophe develops. 

tangential even in tlie absence of tliermal conduction and 
the HBI. However, this tangential velocity bias is also 
lifted for FrMHD ^ 0(1). This has several immediate 
consequences: 

• The boost in thermal conduction restores the heat 
flow from the hot outer layers to the central cool 
core. This is crucial to averting massive overcooling 



Fig. 7. — Time evolution of the Spitzer fraction (solid line) for 
the same run as in Figure 6. Unlike in the HBI case (dashed line), 
the Spitzer fraction, even with cooling, approaches close to / ~ 0.3. 



and heating the cluster cores in a stable fashion. In- 
deed, our runs with radiative cooling showed that a 
cooling catastrophe was averted even in the absence 
of AGN heating. Thus, naive models which assume 
thermal conductivity with Spitzer fraction / ~ 1/3 
may perhaps be reasonably applied to clusters after 
all. Since conduction is regulated by turbulence, it 
will in general be time-dependent. A sudden boost 
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in thermal conduction can modulate a transition 
between a cool core (CC) to a non cool core (NCC) 
cluster (IGuo k OhI 1200^ : we also see some hints 
of this in Fig [7]). Thus, mergers, while generally 
insufficient in the pure hydrody namic case to effec t 
a CC to NCC metamorphosis (jPoole et alJl2008f ). 
may be effective once the boost in thermal conduc- 
tion due to additional turbulence is taken into ac- 
count. If AGN are the main source of turbulence, 
this implies feedback between accretion onto the 
AGN and thermal conduction. Thus, AGN could 
exert another level of self-regulation and thermo- 
static control beyond straight kinetic heating. 

• Our predictions might be testable in the future: 
magnetic topology can be measured by radio po- 
larization measurements (e.g., Pfrommer & Dursi 
2009), while turbulence can be measured by a high 
spectral-resolution calorimeter. Alternatively, our 
predictions can be used to provide indirect con- 
straints on turbulence. For instance, the absence 
of large-scale ordered fields would require a mini- 
mal level of turbulence, and argue against a fully 
relaxed ICM. The resulting turbulent pressure sup- 
port affects mass estimates which assume hydro- 
static equihbrium (L au et al.,i2009,) . and the use of 
clusters for cosmology. 

• The distribution of metals in clusters is known to 
be broader than that of the stars. As shown by 
iSharma et all (|2009allB . metal mixing in a strati- 
fied plasma is much more effective once conduction 
is at play, because of the weaker nature of restor- 
ing buoyancy forces. This allows broad mctallicity 
profiles to be obtained without inversion of the cen- 
tral entropy profile (which is not observed). The 
restoration of thermal conduction by stirring mo- 
tions greatly aids this process. In the limit where 
conduction is so efficient that the cluster becomes 
isothermal (so that 0) or conduction be- 
comes isotropic, fiuid elements have the same tem- 
perature (and hence density) as their surrounding, 
and fiuid elements become neutrally buoyant. This 
makes metal mixing maximally efficient. We note 
that any mechanism which disperses metals over 
wide distances (e.g., turbulent diffusion, Rebusco 
et al. 2008) generally excites turbulent velocities 
sufficient to isotropize the magnetic field. Thus, 
the relatively broad metallicity profiles in clusters 
is consistent with non-negligible stirring. 

Our stirring methods ensure by construction that the 
ensuing turbulence is volume-filling. However, this is by 
no means guaranteed in nature. For instance, galactic 
wakes will be relatively narrow, with a cross-section of or- 
der the galaxy size; similarly, rising AGN bubbles will not 
induce velocity fluctuations over a large solid angle. A 



straightforward way to ensure volume-filling turbulence, 
as is required to stem the HBI, is to excite trapped g- 
modes which are repeatedly refiected and focused within 
a resonance region where uj < co^y^, as discussed in fj2l 
At the same time, oj > ^bv^^ required to overwhelm 
buoyancy forces, an apparently contradictory require- 
ment. Of course, the assumption of a single frequency 
is too simplistic: stirring motions contain many harmon- 
ics and will be scale dependent, generally increasing in 
frequency as turbulence cascades toward smaller scales. 
Thus, large scale g-modes could fall below the Brunt- 
Vaisala frequency, ensuring trapping and amplification, 
but the resulting small scale turbulence could potentially 
still have turnover frequencies exceeding uj^y^. We will 
examine this in future work. 

Although we have focused on the HBI in this paper, 
the same Froude number considerations apply with 
equal force to the MTI unstable outer regions of the 
cluster, where temperature falls with radius and the 
MTI causes field lines to become preferentially radial. 
Here, if Fr^^^ > 0(1), turbulent motions can also 
overwhelm buoyant forces and randomize the field. 
Recently, there has been tantalizing evidence from the 
polarized emission surrounding the magnetic drapes of 
galaxies sweeping up magnetic fields in the Virgo cluster, 
that outside a central region, the magn etic field might 
be p referentially oriented radially (Pfro mmer fc Dursil 
as one might expect from the MTI. If this result 
continues to hold up, this would imply the Froude 
number is below the critical value at these radii (thus 
providing an indirect constraint on turbulent velocities) , 
or additional physical processes, not modeled here, are 
at play. 
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